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A,  ,  " 

I c ;i t  he  results  of  a  Monte  <.  aii>  •»»  udy  of  t  he  most  commonly  discussed 
hods  for  constructing  approximate  confidence  regions  and  confidence  intervals  for 
imcters  estimated  by  nonlir.-ar  least  squares.  The  methods  vie  examine  are  the 
e  variants  of  the  linearization  nvthod,  the  likelihood  method,  and  the  lack-of-Ct 
itod.  T h"  linearization  method  is  the  most  frequently  implemented  method.  It  is 
puiationally  inexpensive  and  produces  easily  understandable  results.  The  likeli- 
i  and  lack-of-Ct  methods  both  are  much  more  expensive  and  more  difficult  to 
rt .  Based  on  our  results,  we  conclude  that  among  the  three  variants  of  the  linear)- 
,i  met  tied,  the  variant  based  solely  on  the  Jacobian  appears  preferable  because  it 
nipler.  le.'S  expensive,  more  numerically  stable,  and  v  least  as  accurate  as  the 
■  t  wo  variants  which  utilize  the  full  Hessian.  !d  our  tests,  however,  all  three  vari¬ 
ed  the  linearization  method  often  produce  gross  underestimates  of  confidence 
>«!'  and  sometimes  produce  significant  underestimates  of  confidence  intervals 
i  the  likelihood  and  lack-of-fit.  methods,  on  the  other  hand,  perform  very  reliably. 

:  tie  datasets  analyzed,  the  Bates  and  Watts  curvature  measures  reliably  predict 
i  the  linearization  method  confidence  regions  will  be  poor,  and  for  the  most  part 
■  s- -  i  -tent  with  our  re  -  ut'  ■  for  the  likelihood  met  had 
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1.  Introduction 


This  paper  presents  the  results  of  an  empirical  study  comparing  several  methods  for  con¬ 
structing  confidence  regions  and  confidence  intervals  about  parameters  estimated  by  nonlinear 
least  squares.  The  methods  compared  are  the  lack-of-fit  method,  the  likelihood  method,  and 
three  variants  of  the  linearization  method. 

The  need  for  confidence  regions  and  intervals  commonly  arises  in  data  fitting  applica¬ 
tions,  where  a  response  variable  y,  observed  with  unknown  error  e,  is  fit  to  m  fixed  predictor 
variables  x,  using  a  function  /(x(;0)  which  can  be  either  linear  or  nonlinear  in  the  p  parame¬ 
ters  0.  The  function  /(x,;0)  is  linear  in  0  if  it  can  be  written 

/(x;;0'  =  x,0  =  ]£)z,;0;,  i’=  1  ,...,n. 

Otherwise,  it  is  nonlinear.  The  methods  analyzed  in  this  study  are  identical  when  /(x,;0)  is 
linear  in  0;  otherwise  they  are  not. 

When  the  error  e,  is  additive,  the  response  variable  can  be  modeled  by 

y,  “  /(x,;0)+e„ 

w  here  0  denotes  the  true  but  unknown  value  of  the  parameters.  The  least  squares  estimator 
of  0  is  the  parameter  value,  denoted  0,  which  minimizes  the  sum  of  the  squares  of  the  residu¬ 
als.  where  the  residuals,  r,(0),  are  estimates  of  the  random  error,  e,, 

',(•)  *  y,-7(x,;0). 


Thus. 


0  =  arg  min  5(0) 


where  5(0)  is  the  residual  sum  of  squares, 

5(«)»  2r,(0)2-  R(0)rR(0) 

i-i 

with  R(0)  denoting  a  column  vector  with  component  r,(0),  and  R(0)r  denoting  the  tran¬ 
spose  of  R(0). 

in  our  study,  we  assume  that  the  model  is  correct  and  that  the  errors  are  normal, 
independent,  identically  distributed  random  variables  with  zero  mean  and  variance  d2,  i.e., 
distributed  as  ,V(0,<r2  I).  Then,  the  least  squares  estimator  0  is  the  maximum  likelihood  esti¬ 
mator  of  the  parameters  0  of  the  p-variate  normal  density  function, 


<> 


A(Y)  =  (2imx-)-"/2  e*~ e ' 

u  Y  is  a  column  vector  with  ilh  component  */,,  and  e  is  a  column  vector  with  ith  com¬ 
ponent  e,. 

Nearly  normally  distributed  errors  are,  in  fact,  encountered  quite  frequently  in  practice. 
This  is  because  measurement  errors  are  often  the  sum.  of  a  number  of  random  errors  from  unk¬ 
nown  sources,  and,  by  the  central  limit  theorem,  the  sum  of  these  errors  is  approximately  nor¬ 
mally  distributed  whatever  the  distribution  of  the  individual  errors  that  make  up  the  sum. 

In  practice,  the  estimated  values  of  the  parameters  0  will  not  equal  the  true  values  0 
because  of  the  random  errors,  e,.  in  the  data.  Since  0  is  a  r -ndom  variable,  howe'er,  it  may 
be  possible  to  indicate  with  some  .specific  probability  I  — a  in  whit  region  about  0  we  might 
reasonably  expect.  0  to  be.  Such  regions  are  known  as  100(1  — ot)r7  confidence  regions.  A  joint 
confidence  region  about  all  of  the  parameters  is  defined  using  a  function 

CRa  :  Y-  a  region  in  Rp 

which  satisfies 

Pr\  e  €  CRJY)  ]  =  1-ct. 


Similarly,  a  confidence  interval  about  an  individual  parameter  0;  is  defined  using  a  function 

C'l,  a  :  Y-  an  interval  in  R 


w Inch  sat isfies 


Pr[  9,  €  Cl,JY)  j  =  l-ot 


The  above  definitions  state  that,  before  the  data  are  sampled,  the  probability  that  the 
confidence  regions  and  confidence  intervals  to  be  constructed  will  contain  the  true  parameter 
values  is  1  — a.  Thus,  if  we  repeatedly  draw  samples  arm  construct  confidence  regions  and 
inter'  ah  about  the  least  squares  estimates  for  each  sample,  then  in  the  long  run  100  (1  —  a)f'c 
of  these  confidence  regions  and  intervals  should  contain  the  true  values.  Methods  that,  for  all 
f  u  net  ion  -  /(x,:0)  and  confidence  levels  1  — a.  are  statistically  guaranteed  asymptotically  to 
contain  the  true  value  100(1  —  a  )rc-  of  the  time  are  called  exact:  all  other  methods  are  called 
approximate. 

Various  methxls  have  been  [imposed  for  calculating  confidence  regions  and  intervals  for 
parameter  estimation  by  nonlinear  least  squares.  These  include  several  variants  of  the  lineari- 
ration  method,  as  well  as  the  likelihood  and  lark-of-fit  methods  [Sec  e.g.  Hard  (107  1).  Gal¬ 
lant  (1976),  Draper  and  Smith  (1081).]  We  review  all  these  methods  briefly  in  Section  2 .  They 
all  arc  equivalent,  and  exact,  for  linear  model*.  For  nonlinear  models,  only  the  lark-of-fit 
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method  for  computing  confidence  regions  is  exact;  the  other  methods  for  computing 
confidence  regions  and  all  the  methods  for  computing  confidence  intervals  are  approximate. 
The  linearization  regions  and  intervals  appear  to  be  the  most  approximate  for  nonlinear 
models,  hut  they  also  are  far  less  expensive  to  compute  than  the  likelihood  or  lack-of-fit 
regions  and  intervals,  and  are  the  predominant  methods  implemented  in  production  software. 
Some  nonlinear  least  squares  packages,  including  NL2SOL  [Dennis,  Gay,  and  Wclsch  (1981)], 
include  three  variants  of  the  linearization  method,  which  differ  only  in  that  the  variance- 
covariance  matrix  of  the  estimated  parameters  is  approximated  in  three  different  ways, 
namely 

va  =  i2  (j(e)rj(§))"\ 

V4  =  s2  H(0)-1, 
or 

vc  =  s'  H(0)- 1  (j(e)TJ(6))  H(S)->, 

where  s-  =  5(0)/(n-p)  is  the  estimated  residual  variance;  J(0)  is  the  Jacobian  of 
/(x,;0),  i-  at  0;  and  H(0)  is  the  Hessian  of  S(0)  at  8. 

Sections  3-6  of  this  paper  describe  and  analyze  a  Monte  Carlo  study  that  compares  all  of 
these  methods  for  computing  confidence  regions  and  intervals  on  20  nonlinear  models.  The 
study  is  used  to  empirically  observe  how  often  the  true  parameter  values  are  contained  in  the 
confidence  regions  and  confidence  intervals  constructed  using  a  given  method.  The  actual 
percent  of  the  nominally  100  (1  —  a)%  confidence  regions  and  intervals  which  are  found  to 
contain  the  true  value  is  known  as  the  observed  coverage.  The  observed  coverage  will  gen¬ 
erally  depend  on  the  method  used  to  construct  the  confidence  regions  and  confidence  inter¬ 
vals.  on  the  nominal  confidence  level,  1-a,  on  the  degree  of  nonlinearity  of  the  function, 
/(x,:0),  and  to  a  small  extent,  on  the  number  of  replications  in  the  simulation.  If  the  experi¬ 
ment  used  to  generate  the  data  is  repeated  a  large  number  of  times  under  the  same  condi¬ 
tions,  and  if  C'Ra  and  Cl  a  are  exact  and  the  model  is  correct,  then  the  observed  coverage  will 
approach  the  nominal  coverage.  When  CR„  and  C7;  a  are  only  approximate,  (he  observed  cov¬ 
erage  will  not  necessarily  approach  the  nominal  coverage,  although  one  would  hope  that  the 
difference  between  the  observed  and  nominal  coverage  for  a  reasonable  approximate  method 
would  be  small  for  most  functions. 

No  similar  study  of  this  magnitude  appears  to  have  been  reported  previously.  The  pro¬ 
perties  of  confidence  regions  and  confidence  intervals  computed  using  the  linearization, 
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likt-lil.u  id,  i’.:!  i  i<  k  i-i'  !ii  im  t  hods  have  bee.  analyzed  by  -evcal  author''.  including  .!<•  nu  rich 
I  1  9-  9)  ! '  ■  ■  r  i  i  ■  •  (J'.ti'a 1 1>;  i  tin  an  and  Mceter  ( 1905),  (i  all  ant  (197b).  I'uncan  I197.H),  and  Hates 
and  W  tils  (  I  9vq)  i.  While  tin-  literature  includes  nu  merous  warnings  regarding  the  possible 
inaccuracy  of  the  approximat  e  method",  it  contains  liitle  empirical  data  '<>  illustrate  the  size 
of  the  discrepancies  between  observed  and  nominal  coverage  that  might  be  expected.  In  those 
studies  which  do  contain  empirical  lata  or.  confidence  re.-'an-  and  intervals,  the  largest 
reported  differences  between  the  observed  and  nominal  coverage  is  only  9rr  for  a  9 5 rc 
confidence  region  computed  using  the  linear/nt  ton  method,  and  is  even  smaller  for  flic  likeli¬ 
hood  method  [Gallant  (1970)].  In  many  practical  applications,  potential  differences  of  9M: 
might  not  be  cause  for  concern.  Kvideure  of  much  larger  differences,  however,  would  indicate 
the  need  for  improved  methods.  Ottr  results  provide  such  evidence. 

Our  Monte  Carlo  study  has  several  purposes.  First,  we  wish  to  determine  whether  the 
observed  coverage  of  the  linearization  method  is  significantly  affected  by  how  the  variance- 
covariance  matrix  is  computed.  Second,  we  wish  to  determine  whether  the  approximate 
confidence  regions  and  confidence  intervals  constructed  using  the  linearization  and  likelihood 
methods,  and  the  approximate  confidence  intervals  constructed  using  the  lark-of-fit  method 
have  observed  coverage  significantly  different  from  nominal.  In  particular,  w want  to  know 
whether  the  frequently  used  linearization  method  is  significantly  better  or  worse  than  the 
more  expensive  likelihood  and  lack-of-fit  methods.  Section  3  describes  how  we  designed  our 
study  to  answer  these  questions.  The  results  are  presented  and  discussed  in  Section  4.  We 
have  also  investigated  how  effective  the  diagnostics  of  Bates  and  Watts  (1980)  are  in  predict¬ 
ing  when  the  confidence  regions  produced  by  the  linearization  and  likelihood  methods  should 
be  reliable:  this  part  of  the  study  is  the  subject  of  Section  5 

Our  study  is  oriented  toward  nonlinear  least  squares  software  developers  who  need 
assiirano'  that  the  methods  they  implement  are  reasonable  for  a  wide  variety  of  problems. 
We  make  only  the  customary  assumptions  that  the  model  is  correct  and  that  the  errors  are 
normally  distributed.  We  do  not  asinine  that  we  can  change  the  representation  of  the  param¬ 
eters.  e.g.,  by  reparameter, 'zing  0  as  log(0).  in  order  *o  reduce  the  difference  between  the 
observed  and  nominal  coverage,  because  reparameterization  is  not  a  technique  that  can  be 
routinely  implemented  by  sofiware  developers  who  have  no  control  ovr  ?  lie  functions 
analyzed.  Headers  interested  in  using  reparameterization  to  improve  their  re-n'ts  are  refered 
to  I .'  at  k  owsky  (  1  9S8 ). 


The  conclusions  we  draw  from  tins  study  are  presented  in  Section  9.  The  first  conclusion 
is  that  among  the  variants  of  the  linearization  method,  the  one  using  Va  is  ilie  besr  choice 
because  it  is  the  cheapest,  and  is  always  at  •. •  *  as  reliable  as  the  other  two  variants  and 
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sometime*  more  reliable.  The  second  conclusion  is  that  even  the  best  linearization  method 
can  be  very  poor;  confidence  regions  with  observed  coverage  as  low  as  12.1%  for  a  nominal 
95%  region,  and  confidence  intervals  with  observed  coverage  as  low  as  75.0%  for  a  nominal 
95%  interval  are  reported.  In  contrast,  for  each  of  the  datasets  tested,  the  confidence  regions 
and  confidence  intervals  constructed  using  the  likelihood  method  and  lack-of-fit  methods  are 
cpiite  close  to  nominal.  Finally,  our  study  indicates  that  the  diagnostics  of  Bates  and  Watts 
(1980)  appear  quite  successful  at  predicting  when  linearization  confidence  regions  will  be  poor. 
Our  recommendations  as  to  how  nonlinear  least  squares  software  should  calculate  confidence 
regions  and  intervals,  in  light  of  these  conclusions,  also  are  given  in  Section  0. 


2.  Background 

This  section  briefly  discusses  methods  for  constructing  confidence  regions  and  confidence 
intervals.  First,  we  give  a  very  quick  survey  of  confidence  regions  and  confidence  intervals  for 
linear  least  squares.  Next,  we  describe  the  two  different  ways  function  nonlinearity  can  affect 
the  solution  locus.  Then,  we  review  the  linearization,  likelihood,  and  lack-of-fit  methods  for 
constructing  confidence  regions  and  confidence  intervals  when  the  model  is  nonlinear.  For  a 
more  complete  discussion,  see  Bard  (1974),  Gallant  (1976),  Draper  and  Smith  (1981),  or 
Donaldson  ( 1985). 

Linear  least  squares 

When  /(x,;9)  is  linear  in  the  parameters  0,  then  /(x,;0)  =  x,  0.  Consequently,  the  Jaco¬ 
bian  of  F(0)  is  X.  an  n  by  p  matrix  with  ith  row  x,.  If  we  assume  that  X  is  of  full  rank,  then 
XrX  is  nonsingular.  and  the  linear  least  squares  estimators  can  be  expressed  in  closed  form 
by 

0  =  (XrX)-;  xty  . 

When  e~.V(0,<r-  I),  a  100  (l-ot)%  confidence  region  about  0  contains  those  values  0  for 
which 

S(I)-5(8)SrpfM.M.„.  (2.1) 

Equation  (2.1)  is  equivalent  to 

(0-0)rXrX(0-0)  £  rpV,.i-. 


(2.2) 


f . , r  ;,|l  ii [i<  r  m  w  i.i.  h  -hows  that  the  shape  of  the  confidence  regions  about  0  ts  ellip¬ 

soidal  f ' ' i  'til  i:nt  ar  n  o  dels 

A  1  ('•>')  i  ’  --  a  j‘V  coniiii' nee  interval  aho.it  6 ,  cotiitnus  i ‘.ease  values  0  fer  which 

s  ’  Vixi,,;1  tn-p.i-a/ .• 


C2.3) 


here  iX~X);i"!  is  the  clement  of  tee  inverse  of  X7X.  The  limits  of  this  confidence 


inter'. d  can  be  shown  to  be  those  values  0;  which 


maximize  (6,-8,)"  subject  to 
5(0)  -  5(0)  =  .1-  (fa.p  -  3-  Fu„. 


>.A) 


p,\- a 


Nonlinearity  and  the  Solution  Locus 

The  solution  locus,  or  estimation  space,  of  /(x,;0 ),  » -  1 . n.  consists  of  all  points  with 

coordinates  expressible  as 

(/(x,;0),/ix  :8) . /<xn:0)j 

whc-e  the  x,.  t  =  I . ri .  are  the  fixed  values  of  the  predictor  variables,  and  0  is  allowed  to 

van  over  all  possible  values  of  the  p  unknown  parameters.  The  solution  locus  is  planar  if 
!  h  r<  e  .1-'  -  a  rep  ar,.  .'ri  z  a*  ion  of  /  (  x , .  8  )  t  ha!  makes  the  function  linear  in  'he  p  parame- 
<  •' i  in  r w  i-e.  t  ho  solut  ie.n  loc us  is  curved. 

\  coordinate  grid  on  the  solution  lorn-  ran  be  formed  by  tricing  the  paths  obtained 
w  |.e  a  e  mh  parameter  i-  individually  allowed  to  vary  while  all  other  parameters  are  held  fixed. 

he  coordinate  grid  i-  curvilinear  w  lie  never  the  function  /( x,  0 )  is  nonlinear  in  one  or  more  of 
it-  pe.r  -rs  ji  i-  lui-ar  only  when  the  'unciion  itself  is  linear. 

ox.-  f  t !;.  -  ]  i.i  i .  >  1 1  l-.-iis  i-  called  ’’intrinsic’’  curvatnr*  [He  de  ( lOOOi:  Hales  and 

V',  c  • '  i  i .  > '  •  ■  i  i.  urv  »•  u:  ••  nl  the  coordinate  grid  is  called  parainef  er-effe  'tv  or  simply 
'par  nii''i  er"  •  ;r'a  ui’  i mi  'A  aits  ( i  rs-0 ) : .  Intrinsic  c  e*--.  .,t  ure  is  mu  alfeeted  by 
reparamet  >Ts/a!  K'n  .ratne*  er--  l”vt«  curvature  is.  linear  function--  have  zero  parameter- 
e  iTee  t  -  curv  at  ure  and  o-ro  it:'':;.'  r  r  v  a*  nr  ••  Noniim  ar  fum"  ems  always  have  nonzero 
paratti' ter-eib  et-.  c  irv  .v  are.  uid  have  color  zero  or  noizero  in*  rui-ie  curvature,  i.e  .  a 

planar  x  -I.’-’'-.!  sol:.'  ;•  a  to.  u  -  el ; 
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Nonlinear  Least  Squares 

When  the  function  is  nonlinear,  the  least  squares  estimators  of  the  parameters  cannot  in 
general  be  expressed  in  closed  form,  and  must  instead  be  computed  by  iterative  techniques. 
Construction  of  exact  confidence  regions  and  confidence  intervals  also  is  much  more  difficult, 
and  so  approximate  methods  are  frequently  used.  The  leading  methods,  linearization,  likeli¬ 
hood,  and  lack-of-fit,  are  described  briefly  below. 

Linearization  methods.  Linearization  methods  for  constructing  confidence  regions 
and  confidence  intervals  assume  that  the  nonlinear  function  can  be  adequately  approximated 
by  an  affine,  or  linear,  approximation  to  the  function  at  the  solution.  That  is,  this  method 
assumes  that  the  solution  locus  is  planar,  and  that  the  coordinate  grid  is  linear  throughout 
the  area  to  be  covered  by  the  confidence  regions  and  confidence  intervals.  Under  this  assump¬ 
tion,  linear  least  squares  theory  tells  us  that  the  confidence  region  about  0  consists  of  those 
values  0  for  which 


(0-0)rV_1(®-®)  s  PFp  n.p  l.a 

while  a  confidence  interval  about  9;,  j—  l . p,  consists  of  those  values  9;  for  which 


where  V  is  the  estimated  variance-covariance  matrix  of  the  parameters,  and  V;}  is  the  (/,/)** 
element  of  V. 

Three  approximations  to  V  are  frequently  used.  These  are 


V3  =  ,-(j(0)rj(0))'\ 

(A) 

V4  =  s-  H(0)_1, 

(B) 

and 

Vc  =  J-H (0)-1  (j(0)rJ(0))  H(0)-‘, 

(C) 

where  J(0)  is  the  Jacobian  of  F ( 0 )  at  0;  H(0)  is  the  Hessian  of  •‘>(0) 
variance.  )2  =  5(0)/n  —  p  . 

at  0;  and  j-  is  the  residual 

Approximation  (A)  is  the  most  common  approximation  to  V. 
mating  F(0)  by  the  affine  approximation  around  0, 

It  is  computed  by  approxi- 

F(0)  ~  F(0)+ J(0)  (0-0)  . 

where  F(0)  is  a  column  vector  with  »(*  component  /(x,;0),  and  then  directly  applying  the 
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intervals  requires  the  solution  of  a  scries  of  nonlincarly  constrained  optimization  problems;  it 
may  be  necessary  to  construct  special  purpose  software  to  solve  these  problems  as  efficiently 
as  possible.)  Performing  hypothesis  tests  using  the  likelihood  or  lack-of-fit  methods  is  compu¬ 
tationally  simple  for  both  confidence  regions  and  intervals,  so  we  recommend  that  one  of  these 
two  methods  be  employed  for  hypothesis  tests  whenever  possible. 

Users  may  prefer  the  likelihood  method  to  the  lack-of-fit  method  even  though  it  is 
approximate  and  the  lack-of-fit  method  is  exact,  because  the  likelihood  method  has  more 
desirable  structural  characteristics  than  the  lack-of-fit  method.  Our  study  provides  no  empir¬ 
ical  evidence  that  the  results  produced  by  the  likelihood  method  are  inferior  to  those  pro¬ 
duced  by  the  lack-of-fit  method.  'Phis  does  not  guarantee  that  similar  results  will  be  obtained 
on  other  datasets,  however.  In  particular,  the  results  of  the  diagnostic  test  proposed  by  Bates 
and  Watts  showed  that  all  our  datasets  have  low  intrinsic  curvature,  which  is  precisely  the 
situation  when  likelihood  methods  are  expected  to  be  very  reliable.  The  additional  dataset  we 
analyzed  with  high  intrinsic  curvature  also  produced  likelihood  method  confidence  region 
observed  coverage  close  to  nominal.  Additional  analysis  is  required  to  determine  whether  the 
hkenhood  method  is  reliable  for  datasets  with  high  intrinsic  curvature,  and  to  determine 
whether  the  Bates  and  Watts  measure  of  intrinsic  curvature  is  a  useful  tool  for  indicating 
when  the  likelihood  method  confidence  regions  are  likely  to  be  unreliable. 

In  addition  to  diagnostics,  it  appears  that  there  is  a  need  for  new  methods  for  estimating 
confidence  regions  that  are  both  reliable  and  easy  to  report.  We  are  especially  interested  in 
investigating  two  methods  that  would  result  in  conservative  elliptical  confidence  regions.  The 
first  method  is  to  find  the  minimal  magnification  of  the  linearization  confidence  region 

that  encloses  the  (9orc)  likelihood  or  lack-of-fit  confidence  region.  This  would  require  the 
solution  of  a  constrained  optimization  problem  with  one  nonlinear  equality  constraint.  The 
second  method  is  to  find  the  smallest  volume  ellipse  that  encloses  the  desired  likelihood  or 
lack-of-fit  confidence  region.  This  would  require  the  solution  of  a  semi-infinite  programming 
poblem.  i  e  an  optimization  problem  with  an  infinite  set  of  constraints. 


7.  Summary 

We  have  presented  the  results  of  a  Monte  Carlo  study  comparing  the  linearization,  likeli¬ 
hood  and  lack-of-fit  methods  for  constructing  confidence  regions  and  confidence  intervals.  Our 
results  indicate  that  the  linearization  method  should  be  constructed  using  the  simplest 
approximation  to  the  variance-covariance  matrix,  (6.1),  as  it  is  simpler,  less  expensive,  more 
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and 

Vc  =  4=  H(0)-»  (j(0)7J(e))  H(0)->  ,  (0.3) 

respectively.  Variant  A  is  simpler  and  less  expensive  because  it  only  requires  the  Jacobian  of 
the  model  function  at  the  solution  and  nol  the  additional  second  order  terms  that  are  also 
required  to  form  the  Hessian.  It  is  more  stable  because  it  can  be  formed  by  inverting  the 
upper  triangular  factor  R  of  the  QR  factorization  of  the  Jacobian  rather  than  by  calculating 
the  inverse  of  the  Hessian;  the  former  calculation  can  be  expected  to  lose  roughly  half  ns  many 
digits  as  the  latter  in  finite  precision  arithmetic. 

The  linearization  method  is  not  always  an  adequate  method  for  approximating 
confidence  regions  and  confidence  intervals  for  the  parameters  of  a  nonlinear  model,  however. 
The  results  presented  in  the  preceding  section  show  just  how  poor  the  linearization  method 
can  be  in  some  cases.  Although  there  are  many  examples  where  the  linearization  method’s 
observed  coverage  differs  from  nominal  by  only  a  very  small  amount,  there  are  also  many 
cases  where  the  observed  coverage  is  far  lower  than  the  nominal.  In  our  tests,  the  best  lineari¬ 
zation  method  variant.  A,  produced  observed  coverages  as  low  as  12. 4^  for  nominal  95% 
confidence  regions  and  75. for  nominal  95^  confidence  intervals. 

Users  will  continue  to  use  the  linearization  method,  however,  because  it  is  readily  avail¬ 
able  in  software  packages  and  provides  a  concise  representation  of  the  information  needed  to 
construct  confidence  regions  and  intervals.  The  erratic  results  obtained  in  our  study  when 
using  the  linearization  method  lead  us  to  conclude  that  users  of  nonlinear  least  squares 
software  must  be  helped  to  cautiously  assess  the  results  they  obtain  using  the  linearization 
method.  The  results  of  the  preceding  section  show  that  the  diagnostic  tools  proposed  by  Bates 
and  Watts  (1980)  are  very  successful  in  indicating  cases  where  the  linearization  method 
confidence  regions  are  likely  to  be  unreliable.  In  these  cases,  more  reliable  methods,  such  as 
the  likelihood  or  lack-of-fit  methods,  are  required  to  produce  accurate  confidence  regions  or 
intervals. 

Our  study  shows  that  the  lack-of-fit  and  likelihood  methods  both  produce  observed  cov¬ 
erages  acceptably  close  to  nominal  in  every  test  case.  Although  the  difficulties  and  expense 
associated  with  using  these  two  methods  to  compute  confidence  regions  make  it  unlikely  that 
they  will  ever  routinely  replace  the  commonly  used  linearization  method  for  this  purpose, 
they  appear  to  be  a  reliable  alternative  that  should  be  considered  when  diagnostics  show  that 
linearization  confidence  regions  are  unreliable.  It  is  not  as  difficult  and  expensive  to  construct 
confidence  intervals  using  the  lack-of-fit  or  likelihood  methods,  and  we  believe  that  producers 
of  nonlinear  least  squares  software  should  consider  this  possibility.  (Constructing  these 
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Cook,  Tsai  and  Wei  (1984)  provide  an  example  which  has  scaled  parameter  effects  curva¬ 
ture  of  934.5  and  scaled  intrinsic  curvature  of  8.4.  Both  the  parameter  effects  curvature  and 
intrinsic  curvature  of  this  dataset  exceed  any  curvature  measure  we  observed  in  the  20 
datasets  in  our  study.  For  this  dataset,  we  computed  observed  confidence  region  coverages  of 
19.0%  and  95.0%  using  the  linearization  method  and  likelihood  methods,  respectively.  While 
the  linearization  method  confidence  region  observed  coverage  is  very  far  from  nominal  as  we 
would  expect  based  on  the  parameter  effects  curvature  of  this  model,  the  likelihood  method 
confidence  region  observed  coverage  is  not.  We  cannot  conclude  anything  from  this  one  obser¬ 
vation.  It  is  clear,  however,  that  additional  analysis  of  datasets  with  high  intrinsic  curvature 
would  be  useful  to  further  assess  the  effect  of  a  non-planar  solution  locus  on  the  likelihood 
method. 


6.  Conclusions 

Based  on  our  computational  study,  we  can  draw  conclusions  about  :  i)  the  comparison 
between  the  three  variants  of  the  linearization  method;  ii)  the  reliability  of  linearization 
methods  for  calculating  confidence  regions  and  confidence  intervals;  and  iii)  the  reliability  of 
the  likelihood  and  lack-of-fit  methods  for  calculating  confidence  regions  and  confidence  inter¬ 
vals. 

When  using  the  linearization  method  to  construct  confidence  regions  and  intervals,  our 
Monte  Carlo  study  has  shown  no  clearcut  difference  in  the  observed  coverage  of  one  variant  as 
compared  to  another.  In  our  tests,  the  only  statistically  significant  difference  among  the 
results  produced  by  the  three  linearization  variants  was  in  constructing  confidence  intervals 
with  finite  difference  Jacobians  and  Hessians;  here  variant  A  was  superior  to  variants  B  and  C. 
We  found  no  empirical  evidence  that  one  should  prefer  variants  B  or  C,  even  though  they  may 
be  appealing  from  a  theoretical  point  of  view.  Therefore  we  conclude  that  variant  A  of  the 
linearization  method,  which  is  computed  using 

V,  =  t2  (j(0)rJ(9))"1  ,  (6.1) 

is  the  best  variant  to  use  for  constructing  both  confidence  regions  and  confident-'  intervals, 
because  it  is  simpler,  less  expensive,  and  more  numerically  stable  to  compute  than  variants  B 
or  C,  which  use 


V*  -  ^Hfl)"1 


(6.2) 
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locus  is  planar  is  valid  over  the  region  of  interest,  and  therefore  the  likelihood  method 
confidence  region  should  be  adequate. 

In  Figure  5  we  plot  the  20  confidence  region  observed  coverages  obtained  using  lineariza¬ 
tion  method  variant  A  with  analytic  derivatives  (derivative  configuration  DC3)  and 
e~7V( 0,(ti  a)2  i)  against  the  Bates  and  Watts  relative  measure  of  parameter  effects  curva¬ 
ture.  Likewise,  in  figure  6  we  plot  the  corresponding  20  likelihood  method  confidence  region 
observed  coverages  against  the  Bates  and  Watts  relative  measure  of  intrinsic  curvature.  The 
relative  curvature  measures  were  computed  at  the  true  parameter  values  using  the  true  vari¬ 
ance  of  the  errors.  In  these  plots,  we  have  scaled  the  measures  of  parameter  effects  curvature 
and  intrinsic  curvature  by  dividing  the  measure  by  the  appropriate  critical  value.  Thus,  in 
both  of  these  plots,  a  scaled  curvature  measure  less  than  l  indicates  the  relative  measure  was 
less  than  the  critical  value,  while  a  value  greater  than  1  indicates  the  curvature  exceeded  the 
critical  value. 

It  is  clear  from  figure  5  that  the  Bates  and  Watts  parameter  effects  curvature  measure  is 
strongly  correlated  with  the  observed  coverage  obtained  using  the  linearization  method.  In 
fact,  for  our  data,  as  the  parameter  effects  curvature  increases,  the  observed  coverage  for  the 
linearization  method  confidence  regions  decreases  nearly  monotonically  and  linearly  as  the 
logarithm  of  the  scaled  parameter  effects  curvature.  Furthermore,  in  all  datasets  where  the 
parameter  effects  curvature  is  less  than  the  critical  value,  the  observed  confidence  region  is 
very  close  to  nominal,  while  in  all  cases  where  the  parameter  effects  curvature  is  greater  than 
ten  times  the  critical  value,  the  observed  coverage  is  unsatisfactorily  low.  Datasets  with 
parameter  effects  curvature  between  one  and  ten  times  the  critical  value  had  observed 
confidence  region  coverage  between  83.2^  and  91.A?S.  From  these  results,  it  appears  that  the 
Bates  and  Watts  parameter  effects  curvature  is  a  reliable,  if  perhaps  stringent,  indicator  of 
when  the  linearization  method  will  produce  reliable  confidence  regions. 

Figure  6  shows  that  all  but  one  of  the  20  datasets  tested  in  this  study  have  intrinsic  cur¬ 
vature  which  is  less  than  the  critical  value,  which  means  that  each  of  these  datasets  is  nearly 
planar.  For  nearly  planar  datasets  we  expected  good  observed  coverage  from  the  likelihood 
method,  and,  as  figure  6  shows,  that  is  what  we  got.  Since  none  of  our  datasets  have  high 
intrinsic  curvature,  however,  we  do  not  know  how  the  likelihood  method  will  perform  when 
the  solution  locus  is  not  nearly  planar.  We  cannot  assume  that  the  accurate  results  obtained 
in  our  study  using  the  likelihood  method  will  necessarily  carry  over  to  datasets  with  large 
intrinsic  curvature. 


17 


about  this  linear  combination  was  approximately  equal  to  that  of  the  linearization  method 
confidence  region  observed  coverage. 

The  use  of  finite  differences  to  approximate  both  the  Jacobian  and  the  Hessian  appears  to 
significantly  degrade  the  confidence  interval  results  for  linearization  variants  B  and  C.  Figure 
3  shows  that,  while  there  is  no  striking  difference  in  the  results  obtained  using  the  three  vari¬ 
ants  of  the  linearization  method  with  derivative  configurations  DC2  and  DC3,  variants  B  and 
C  degrade  significantly  more  than  variant  A  when  using  DC1,  i.e.,  finite  difference  Jacobian 
and  Hessian.  A  two-sided  paired-sample  t-test  was  used  to  determine  whether,  for  a  given 
derivative  configuration,  the  observed  coverages  obtained  using  the  different  linearization 
method  variants  are  statistically  different  at  the  5%  significance  level.  The  results  indicate 
that  when  derivative  configuration  DC2  and  DC3  are  used,  the  differences  in  the  results 
obtained  using  variants  A,  B,  and  C  are  seldom  statistically  significant  at  the  5 9c  level,  but 
that  when  the  Jacobian  and  Hessian  are  approximated  using  finite  differences  (derivative 
configuration  DCl)  then  the  differences  in  results  are  often  significant. 

Comparing  Figures  3  and  4  shows  that  as  the  variance  of  the  errors  is  increased,  the 
differences  between  observed  and  nominal  coverage  also  are  increased,  as  was  the  case  for  the 
confidence  region  results.  However,  this  increase  is  not  as  pronounced  for  confidence  intervals 
as  for  confidence  regions.  The  results  at  confidence  levels  0.50,  0.75,  0.95,  and  0.99  also 
showed  that  as  the  nominal  confidence  level  approaches  100%,  the  spread  between  observed 
and  nominal  coverages  obtained  using  the  linearization  method  is  increased. 


5.  Diagnostic  tools 

The  preceding  section  demonstrates  a  pressing  need  for  diagnostics  to  warn  users  when 
the  commonly  used  linearization  method  confidence  region  will  not  have  adequate  coverage. 
In  addition,  it  would  be  useful  to  have  a  warning  to  indicate  when  the  approximate  likelihood 
method  may  be  inadequate.  Bates  and  Watts  (1980)  have  proposed  measures  of  nonlinearity 
that  provide  such  diagnostics. 

According  to  Bates  and  Watts,  when  their  relative  measure  of  parameter  effects  curva¬ 
ture  is  small  compared  to  the  critical  value  {Fp  „_p  0  os)~ U2,  then  the  linear  coordinate  grid 
assumption  is  valid  over  the  region  of  interest,  and  therefore  the  linearization  method 
confidence  region  should  be  adequate.  Similarly,  when  their  relative  measure  of  intrinsic  cur¬ 
vature  is  small  compared  to  the  same  critical  value,  then  the  assumption  that  the  solution 
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Confidence  Intervals 

Results.  Figures  3  and  4  provide  information  for  confidence  intervals  which  is  analo¬ 
gous  to  that  shown  in  figures  1  and  2  for  confidence  regions.  The  observed  coverages  plotted 
are  the  smallest  of  the  p  confidence  interval  coverages  obtained  for  each  dataset.  Figure  3 
displays  the  observed  confidence  interval  results  for  nominally  95%  confidence  levels,  when 
e~ /V(0,(r2 1);  figure  4  shows  the  results  when  e~  n(o,(t\  ct)2  i),  excluding  linearization 
method  variants  B  and  C  as  was  done  for  the  linearization  method  confidence  regions. 

Observations.  Figure  3  shows  that  for  confidence  intervals,  the  best  results  are 
obtained  using  the  lack-of-fit  and  likelihood  methods,  and  the  worst  results  are  obtained  using 
the  linearization  method,  as  was  the  case  for  confidence  regions.  The  lack-of-fit  and  likeli¬ 
hood  methods  produce  confidence  intervals  which  seldom  vary  from  nominal  by  an  amount 
that  is  significant  at  the  5%  level,  and  never  are  less  than  nominal  by  more  than  5.0%. 
Again,  use  of  finite  difference  Jacobians  does  not  appear  to  affect  the  results  for  these  two 
methods. 

The  three  variants  of  the  linearization  method,  on  the  other  hand,  frequently  produce  far 
less  reliable  confidence  intervals  than  the  lack-of-fit  and  likelihood  methods.  Disturbing 
differences  between  observed  and  nominal  coverages  occur  when  each  of  the  variants  of  the 
linearization  method  is  used  to  construct  confidence  intervals.  The  observed  coverage  for  a 
nominally  95%  confidence  interval  is  as  low  as  75.0%,  44.0%,  and  10.8%  for  variants  A,  B, 
and  C.  respectively.  For  most  of  the  datasets  tested  in  our  study,  however,  the  span  between 
observed  and  nominal  coverage  produced  by  the  three  variants  of  the  linearization  method  is 
considerably  less  for  confidence  intervals  than  for  linearization  method  confidence  regions  con¬ 
structed  about  the  parameters  of  the  same  dataset.  This  is  especially  true  when  derivative 
configurations  DC‘2  and  DC3  are  used. 

One  reason  why  linearization  method  confidence  intervals  have  better  coverage  than 
linearization  method  confidence  regions  is  that,  when  the  parameter  estimates  are  correlated 
with  each  other,  a  number  of  points  may  be  included  in  the  linearization  method  confidence 
intervals  but  not  in  the  confidence  regions.  Note,  however,  that  if  a  confidence  interval  was 
computed  for  the  linear  combination  of  the  parameters  given  by  the  eigenvector  correspond¬ 
ing  to  the  minor  axis  of  the  linearization  method  confidence  region  ellipsoid,  then  the  lineari¬ 
zation  method  confidence  interval  observed  coverage  should  approximately  equal  that  of  the 
linearization  method  confidence  region.  In  our  Monte  Carlo  study,  we  actually  computed  the 
linearization  method  confidence  iaterval  observed  coverage  for  this  linear  combination  of  the 
parameters.  In  every  case,  the  observed  coverage  we  obtained  for  the  confidence  interval 
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Observations.  Figures  1  and  2  show  that  the  !ack-of-fit  and  likelihood  method 
confidence  regions  arc  quite  reliable,  and  that  the  results  are  not  affected  by  use  of  finite 
difference  derivatives.  In  all  our  tests,  they  produced  observed  coverages  which  seldom  vary 
from  nominal  by  an  amount  that  is  significant  at  the  5%  level.  In  fact,  for  these  datasets, 
there  is  only  one  instance  (dataset  3AAA,  i  ~~  N{o,{r |  d)2l))  where  the  difference  between 
the  nominal  and  observed  coverages  produced  using  these  two  methods  is  greater  than  o%, 
and  in  this  instance,  the  observed  coverage  is  greater  than  nominal,  not  less. 

The  three  variants  of  the  linearization  method,  on  the  other  hand,  frequently  produced 
far  less  reliable  confidence  regions,  although,  as  discussed  below,  the  results  still  do  not 
appear  to  be  affected  by  the  use  of  finite-difference  derivatives.  The  difference  between  the 
nominal  and  observed  coverages  obtained  using  the  linearization  methods  often  are  consider¬ 
ably  more  than  20 %,  which  is  a  difference  that  many  if  not  most  users  would  find  unaccept¬ 
able. 

By  comparing  Figure  1  to  Figure  2,  it  is  apparent  that  increasing  the  variance  of  the 
errors  does,  in  fact,  increase  the  differences  between  observed  and  nominal  coverage  for  all 
methods.  Our  tests  at  confidence  levels  0.50,  0.75,  and  0.99,  which  are  not  reported  in  detail 
here,  also  showed  that  the  spread  between  the  observed  and  nominal  coverage  obtained  using 
the  linearization  method  increases  as  the  nominal  confidence  level  is  increased. 

The  large  differences  for  some  datasets  between  the  observed  coverage  of  confidence 
regions  constructed  using  the  likelihood  method  and  those  obtained  using  the  linearization 
method  may  be  explained  by  the  difference  in  the  shape  of  the  two  regions.  The  likelihood 
method  confidence  region  corresponds  to  the  boundary  and  interior  of  a  contour  of  the  sum  of 
squares  surface,  i.e.,  a  contour  of  constant  likelihood,  whereas  the  linearization  method 
confidence  regions  are  always  ellipsoidal.  We  plotted  these  contours  for  various  datasets,  and 
the  difference  sometimes  were  very  large.  Examples  for  datasets  3AAA  and  N  \AG  are  given 
in  Donaldson  ( 1985). 

Figure  1  also  indicates  that  the  observed  coverage  obtained  using  variants  A,  B,  and  C  of 
the  linearization  method  are  nearly  identical.  The  results  of  two-sided  paired-sample  t-tests 
indicate  that  there  is  no  statistically  significant  differences  at  the  5 °c  level  between  the 
observed  coverages  obtained  using  any  of  the  variants  of  the  linearization  method  with  any  of 
the  derivative  configurations.  The  same  results  were  obtained  for  our  tests  at  the  0.50.  0.75, 
and  0  99  confidence  levels. 
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4.  Results  and  Observations 

This  section  presents  the  results  of  our  Monte  Carlo  study  of  the  lack-of-fit  method,  the 
likelihood  method,  and  the  three  variants  of  the  linearization  method.  The  section  is  divided 
into  a  discussion  of  confidence  regions  and  confidence  intervals.  For  each,  wc  also  make  a 
number  of  observations  about  the  results.  The  conclusions  we  draw  from  our  analysis  are  dis¬ 
cussed  in  the  next  chapter. 

The  material  in  this  chapter  includes  a  number  of  figures.  These  are  printed  at  the  end 
of  the  paper. 

Confidence  Regions 

Results.  The  results  for  nominally  95%  confidence  regions  constructed  using  each  of  the 
methods  analyzed  in  this  study  with  e~/V(0,<x2  I)  are  graphically  displayed  in  Figure  1.  For 
each  dataset,  the  observed  coverage  is  plotted  against  the  method  and  derivative 
configuration  used  to  obtain  it. 

The  three  derivative  configurations  are  labeled  DC1,  DC2,  and  DC3  in  these  and  the  fol¬ 
lowing  figures  and  tables,  as  well  as  in  Appendix  B,  Here  DCl  denotes  use  of  finite  difference 
approximations  for  both  the  Jacobian  and  the  Hessian,  DC2  denotes  use  of  analytic  Jacobian 
and  finite  difference  Hessian,  and  DC3  denotes  use  of  analytic  Jacobian  and  Hessian.  Since 
the  computations  required  to  calculate  the  lack-of-fit  method  results  and  the  likelihood 
method  results  using  derivative  configurations  DC2  and  DC3  are  exactly  the  same,  these 
results  are  displayed  together. 

Figure  2  shows  the  analogous  results  for  e~  /v(o,(si  &)2  i).  As  noted  in  Section  3,  vari¬ 
ants  B  and  C  of  the  linearization  method  are  excluded  from  the  analysis  displayed  in  Figure  2 
because  computational  difficulties  were  encountered  for  these  variants  when  the  variance  of 
the  errors  was  increased. 

A  conservative  95%  confidence  interval  about  the  nominal  confidence  level  is  indicated 
on  each  plot  by  a  pair  of  horizontal  lines  which  represent  the  values  100  ( 1  -  a)±  1. 1,  where 
1. 1  is  two  times  the  maximum  standard  deviation  of  the  observed  coverage  at  any  coverage 
level.  This  confidence  interval  provides  a  quick  means  of  determining  whether  any  of  the 
observed  coverages  for  each  method  are  significantly  different  from  the  nominal  confidence 
level  at  the  5%  level.  When  the  method  used  to  construct  the  confidence  regions  and 
confidence  intervals  is  exact,  we  expect  that  the  observed  coverage  for  9V<  of  all  possible 
datasets  will  lie  within  this  confidence  interval. 
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datasets  arc  often  troublesome,  they  are  mostly  real  world  problems  that  have  not  been  made 
artificially  difficult. 

Each  dataset  was  analyzed  twice  to  allow  us  to  examine  the  effect  of  increasing  the  stan¬ 
dard  deviation  of  the  errors.  In  the  first  analysis,  e  ~  /V(0,c2I);  in  the  second  analysis, 
i  ~  N{0,(r\&yi),  where  ^  is  approximately  the  largest  number  s  10  for  which  every  reali¬ 
zation  of  the  data  could  be  successfully  analyzed.  The  methods  analyzed  in  the  second 
analysis  were  the  same  as  in  the  first  except  that  variants  B  and  C  of  the  linearization  method 
were  excluded  from  the  second  analysis  because,  when  vj>1.0,  we  were  frequently  unable  to 
compute  the  required  test  statistics  using  these  two  variants. 

Computation  of  the  linearization  method  and  the  lack-of-fit  method  requires  that  certain 
derivatives  be  available.  The  Jacobian  of  F(8)  is  used  by  both  the  linearization  and  lack-of-fit 
methods.  Variants  B  and  C  of  the  linearization  method  use  the  Hessian  of  5(8)  as  well.  In 
practice,  analytic  derivatives  often  are  not  available.  Therefore,  in  our  study  each  method 
was  implemented  and  analyzed  using  three  different  derivative  configurations.  These 
configurations  are  (1)  the  Jacobian  and  Hessian  both  approximated  by  finite-differences,  (2) 
the  Jacobian  computed  analytically  and  the  Hessian  computed  by  finite-differences,  and  (3) 
both  the  Jacobian  and  the  Hessian  computed  analytically.  For  derivative  configurations  (1) 
and  (2),  the  variance-covariance  matrix  needed  by  the  linearization  method  was  returned 
directly  from  NL2SOL  (Dennis,  Gay  and  YVelsch  (1981)].  For  configuration  (3),  it  was  con¬ 
structed  outside  of  NL2SOL.  For  details  on  the  formulas  used  to  compute  the  finite-difference 
derivative  approximations,  see  Donaldson  (1985). 

We  ran  our  Monte  Carlo  study  in  single  precision  on  a  60  bit  word  length  computer.  All 
subroutines  extracted  from  other  sources  were  used  without  modification  except  for  NL2SOL, 
which  was  changed  for  this  study  in  two  important  ways.  First  we  disabled  the  two  tests 
within  NL2SOL  used  to  detect  near  singularity.  Second,  we  used  the  STARPAC  [Donaldson 
and  Trvou  (1983)]  front  end  to  NL2SOL.  With  this  front  end.  the  finite  difference  approxima¬ 
tion  to  the  Jacobian  is  computed  with  the  optimal  derivative  step  sizes  selected  using  the 
algorithm  developed  by  Schnabel  (1981),  thus  maximizing  the  number  of  correct  digits  in  each 
element  of  the  finite  difference  Jacobian. 


true  parameter  values.  When  oj> a,  the  true  value  did  not  lie  in  the  100  (1  -  a)%  confidence 
region  or  confidence  interval;  when  it  did.  The  values  1  —  to  were  obtained  using  the 

hypothesis  tests  corresponding  to  the  formulas  for  confidence  regions  and  intervals  given  in 
Section  2,  and  the  appropriate  cumulative  distribution  functions;  the  procedures  arc  described 
in  detail  in  Donaldson  (1985).  The  cumulative  distribution  functions  were  obtained  from  the 
STARPAC  subprogram  library  (Donaldson  and  Tryon  (1983)). 

The  observed  coverage,  7a,  for  the  particular  nominal  confidence  level,  method  and  sys¬ 
tem  under  analysis  is  the  percentage  of  the  total  number  of  realizations  of  the  data,  N,  for 
which  a)  ^  a.  When  N  is  large,  the  standard  deviation  of  /„  can  be  approximated  using  the 
normal  approximation  to  the  binomial  distribution.  In  this  study  we  used  Nm  500,  so  the  max¬ 
imum  standard  deviation  of  the  observed  coverage  at  any  coverage  level  is  approximately 
2.2*0. 

Note  that  substituting  a  new  realization  of  the  data  for  one  which  could  not  be  com¬ 
pletely  analyzed  because  either  (a)  the  nonlinear  least  squares  algorithm  did  not  converge,  or 
(b)  the  test  statistics  could  not  be  computed  for  every  method  being  analyzed,  is  a  form  of 
censoring  which  will  bias  the  observed  coverages  obtained.  In  our  analysis,  we  adjusted  the 
value  of  <r  for  each  dataset  so  that  every  realization  could  be  completely  analyzed,  and  there¬ 
fore  the  results  reported  in  this  paper  are  not  derived  from  censored  data. 

We  computed  the  observed  coverage  for  four  nominal  confidence  levels,  0.50,  0.75,  0.95, 
and  0.99.  In  this  paper  we  only  include  our  data  for  the  level  0.95,  although  we  comment 
briefly  in  Section  4  on  our  results  at  the  other  levels.  Data  for  the  full  study  are  given  in 
Donaldson  ( 1985). 

The  references  for  the  datasets  used  in  our  Monte  Carlo  study  are  given  in  Appendix  A 
and  described  in  detail  in  Donaldson  (1985).  With  only  two  exceptions,  the  functions  and 
data  which  comprise  our  datasets  have  been  taken  from  Ratkowsky  (1983),  Himmelblau 
( 1970).  Guttman  and  Meeter  ( 1965),  and  Duncan  (1978).  The  standard  deviation  of  the  errors 
of  some  of  the  datasets  has  been  adjusted  in  order  to  allow  us  to  successfully  aualyze  each 
realization  of  the  data  for  each  dataset.  The  two  datasets  not  taken  from  the  published 
literature  are  identified  as  8ACA  and  9AAG.  Dataset  8ACA  was  created  especially  for  this 
study  by  generalizing  fuuction  3  to  a  larger  number  of  parameters.  Dataset  9AAG  involves  a 
microwave  absorption  line  function  taken  from  a  consulting  session  at  the  National  Bureau  of 
Standards  in  Boulder,  Colorado. 

The  number  of  parameters  in  the  20  datasets  analyzed  range  from  2  to  8  and  the  ratio  of 
the  number  of  parameters  to  the  number  of  observations  range  from  2/12  to  3/5.  While  these 
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flit  k=  1 1  ,j  +  1  otherwise  it  is  approximate.  [See  Halperin  (1983). j 

The  lack-of-fit  method  is  even  more  expensive  to  use  than  the  likelihood  method,  and,  as 
is  the  case  for  the  likelihood  method,  the  information  needed  to  construct  the  confidence 
regions  cannot  be  succinctly  summarized  for  publication.  Also,  the  confidence  regions  and 
confidence  intervals  constructed  using  the  lack-of-fit  method  are  guaranteed  to  contain  every 
minimum,  maximum,  and/or  saddle  point  of  the  likelihood  surface.  This  makes  the  lack-of-fit 
method  structurally  undesirable. 


3.  The  Monte  Carlo  Study 

This  section  briefly  describes  how  our  Monte  Carlo  study  was  constructed.  Full  details 
are  provided  by  Donaldson  (1985). 

The  Monte  Carlo  method  uses  the  computer  to  simulate  the  results  of  repeating  an 
experiment  many  times  in  order  to  obtain  a  large  sample  from  which  the  statistical  properties 
of  a  system  can  be  examined.  For  each  simulation,  we  first  generated  the  errors  and  response 
variables.  The  errors,  e,  were  produced  using  the  Marsaglia  and  Tsang  pseudo-normal  ran¬ 
dom  number  algorithm  (1984)  as  implemented  by  James  Blue  and  David  Kahanar  of  the 
National  Bureau  of  Standards  Scientific  Computing  Division.  The  response  variable,  Y,  was 
then  constructed  so  that  its  iM  component  is 

y,  ■  /(x,;0)+e,  . 

Then  the  least  squares  estimate,  0,  was  calculated  using  NL2SOL.  an  unconstrained  quasi- 
Newton  code  for  nonlinear  least  squares  [Dennis,  Gay,  and  Welsch  (1981)].  Starting  values  for 
NL2SOL  were  set  to  the  true  values  of  the  parameters,  0,  and  the  stopping  criteria  for  the 
convergence  tests  based  on  the  relative  change  in  the  parameters  and  in  the  sum  of  squares 
both  were  set  to  10-s. 

Finally,  for  each  confidence  region  or  interval  method  and  each  derivative  configuration 
being  analyzed,  we  recorded  whether  the  true  values  of  the  parameters  were  contained  within 
the  confidence  regions  and  confidence  intervals  for  this  realization  of  the  data.  Determining 
whether  the  true  parameter  values  lay  within  the  confidence  regions  and  confidence  intervals 
about  the  least  squares  estimates  fortunately  did  not  require  that  we  construct  the  full 
confidence  regions  and  confidence  intervals  for  each  confidence  level  and  method.  Instead,  we 
simply  calculated  the  smallest  confidence  level,  l-u»,  such  that  a  1 00  (l  —  tu )<~c  confidence 
region  or  confidence  interval  constructed  using  the  method  being  analyzed  will  contain  the 
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intervals  and  confidence  regions  for  subsets  of  the  parameters.  The  lack-of-fit  method  is  based 
on  the  fact  that  the  quadratic  forms 

Ql(e)  =  Mmilsiii 

d2 

and 

g,( .  SI«lrfcFl»})R[e| 

d-2 


where 

P(0)  =  j(e)(j(e)rJ(e))_,j(0)^ 

are  independent  chi-square  random  variables  with  p  and  n-p  degrees  of  freedom  respec¬ 
tively.  Therefore, 

Qi(e)/p 

<?2(0)/(n-p) 

is  distributed  as  Fp  n^p  l_a,  so  an  exact  100  (1  —  a)%  confidence  region  consists  of  all  values  0 
such  that 

— RliilPilMli  f 

R(0)r(l-P(0))R(0)  n~P  p'n~fl~a 

Note  that  the  lack-of-fit  method  does  not  require  that  the  least  squares  solution  be  found 
prior  to  constructing  the  confidence  region. 

Similarly,  a  lack-of-fit  method  confidence  interval  for  the  jth  parameter  consists  of  those 
values  0;  for  which  there  exists  values  of  0t,  *  =  l,j+  l,...,p,  such  that  for  these  p 

parameter  values,  0, 

,y(n-P) 

where  i*  the  residual  sum  of  squares  obtained  when  R(0)  is  linearly  fit  to  all  the 

columns  of  J(0)  excluding  the  jth,  and  S’t(0J( fj)  is  the  residual  sum  of  squares  obtained  when 
R(0)  is  linearly  fit  to  J(0).  This  interval  is  exact  if  /(x,;0)  is  linear  in 
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7,cro,  however,  we  expect  that  the  linearization  methods  could  sometimes  produce  observed 
coverages  very  far  from  the  expected  nominal  coverage.  The  results  of  our  Monte  Carlo  study 
show  this  to  be  true. 

Likelihood  method.  The  likelihood  method  is  another  approximate  method  for  pro¬ 
ducing  confidence  regions  and  confidence  intervals.  The  likelihood  m*  v,od  confidence  region 
about  8  consists  of  those  values  8  for  which 

£(«)-*(«}  s  *-PF,  . 

This  is  analogous  to  equation  (2.1)  for  confidence  regions  for  the  parameters  of  a  linear  func¬ 
tion,  although  when  /(x,;8)  is  nonlinear  in  the  parameters  the  resulting  confidence  region  is 
no  longer  ellipsoidal.  The  likelihood  method  confidence  interval  about  is  the  interval 
bounded  by  the  points  which 

maximize  (8;  — §;)2  subject  to 

s(e)-s(S)ts  f,, 

This  confidence  interval  is  the  projection  onto  the  appropriate  parameter  axis  of  the  above 
region,  and  is  analogous  to  equation  (2.4)  for  confidence  intervals  in  the  case  of  linear  least 
squares. 

When  the  solution  locus  is  planar,  the  confidence  regions  (but  not  the  confidence  inter¬ 
vals)  constructed  using  the  likelihood  method  are  equivalent  to  the  lack-of-fit  confidence 
regions,  and  therefore  are  exact.  In  addition,  likelihood  method  confidence  regions  and  inter¬ 
vals  have  the  desirable  property  that  they  are  constructed  from  contours  of  constant  likeli¬ 
hood,  and  that  the  regions  and  intervals  are  not  affected  by  reparameterization  of  the  func¬ 
tion  /(x,:8).  Thus  we  might  expect  the  likelihood  method  to  produce  confidence  regions  and 
confidence  intervals  with  observed  coverage  closer  to  nominal  than  those  produced  using  the 
linearization  methods.  However,  the  likelihood  method  has  several  practical  disadvantages. 
Both  the  confidence  regions  and  confidence  intervals  produced  using  the  likelihood  method 
can  be  disjoint  and  unbounded  because  the  contours  of  a  nonlinear  function  can  be  disjoint 
and  unbounded.  The  method  also  is  very  expensive  to  use,  and,  when  the  data  arrays  are 
large,  it  can  be  awkward  to  publish  the  information  necessary  to  reconstruct  the  confidence 
region  because  this  information  is  not  succinctly  summarized  as  it  is  in  the  case  of  the  lineari¬ 
zation  method. 

Lack-of-fit  method.  The  lack-of-fit  method  can  be  used  to  produce  exact  joint 
confidence  regions  for  all  p  of  the  parameters,  and  to  produce  approximate  confidence 


'linear  least  squares  theory. 


Approximation  (B)  can  be  obtained  using  maximum  likelihood  theory.  For  large  samples, 
maximum  likelihood  estimators  are  asymptotically  distributed  as  the  p-variate  normal  with 
variances  and  covariances  given  by  V  where 


V-i 


P  fa2logL(Y)l 

ae;ae*  J 


!t  is  stringhtforeward  to  show  that  V-1  approaches  V6  1  as  »-». 

Approximation  (C)  can  be  obtained  from  sensitivity  analysis.  If  the  observations  Y  are 
changed  to  Y  +  t ,  then,  to  within  0(e)  terms,  0  will  be  changed  to 

0(e)  -  0-H(0)-'  J(0)Te. 


Solving 

V  =  coe(0(e))  -  E  ((  6(e)-0  )(  0(e)-0  )*) 


yields  V  =  V£. 

When  certain  regularity  conditions  are  met  [Jennrich  (1959)],  each  of  these  approxima¬ 
tions  to  V  asymptotically  will  approach  the  true  variance-covariance  matrix  of  the  model. 
"■  >!'■  also  that  these  approximations  differ  only  when 


a?/(x,;@) 


s  nonzero.  In  particular,  for  linear  functions,  each  of  these  representations  of  V  is  equal  to 

»2  (j(0)rJ(0))_1  =  j2  (XrX)_1  . 


For  nonlinear  functions,  V6  is  said  to  utilize  observed  information  and  Va  is  said  to  utilize 
expected  information. 


Linearization  methods  have  the  advantage  that  their  resulting  confidence  regions  and 
intervals  are  simple  and  inexpensive  to  construct,  and  that  they  produce  bounded,  convex 
confkl  ence  regions.  In  addition,  the  information  needed  to  construct  confidence  regions  and 
intervals  using  this  method  can  be  parsimoniously  summarized  by  the  p  by  p  matrix  V.  and  is 
well  understood  by  users  familiar  with  linear  least  squares.  Because  the  linearization  methods 
assume  that  both  the  intrinsic  curvature  and  the  parameter-effects  curvature  of  /( x , : 0 )  are 
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numerically  stable,  and  at  least  as  accurate  as  the  other  two  linearization  variants,  which  are 
constructed  using  (0.2)  and  (6.3).  We  have  also  given  considerable  evidence  that  confidence 
regions,  and  to  some  extent  confidence  intervals,  constructed  using  the  linearization  method 
can  be  essentially  meaningless. 

Our  study  shows  that  the  likelihood  and  lack-of-fit  methods,  on  the  other  hand,  pro¬ 
duced  consistently  good  results  for  the  datasets  tested.  However,  because  the  likelihood 
method  is  approximate  it  is  not  clear  that  the  good  results  we  obtained  with  it  will  necessarily 
be  characteristic  of  all  datasets.  Also,  because  of  the  undesirable  structural  characteristics  of 
the  lack-of-fit  method,  it  is  unlikely  to  be  used  routinely,  although  in  cases  where  accuracy  is 
of  extreme  importance,  it  may  be  a  useful  tool  to  have. 

Because  of  the  uncertainty  associated  with  the  linearization  and  likelihood  methods,  we 
also  have  briefly  examined  how  the  Bates  and  Watts  curvature  measures  relate  to  the 
confidence  region  observed  coverages  we  obtained  in  this  study.  Our  results  show  that  the 
Bates  and  Watts  parameter  effects  curvature  appears  to  provide  excellent  indication  of  when 
the  linearization  method  may  produce  less  than  satisfactory  results.  Our  results  are  not  as 
conclusive,  however,  about  the  relation  between  intrinsic  curvature  and  likelihood  method 
coverage  since  the  solution  locus  for  all  of  our  datasets  were  nearly  planar. 
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Appendix 


Dataset  Id.  p/n  Reference 

1  2AAA  2/12  Guttman  and  Meeter  (1965)  ;  model  ,  page  628 

2  3AAA  2/12  Guttman  and  Meeter  (1965)  ;  model  ri3  ,  page  628 

3  4AAA  2/21  Duncan  (1978)  ;  model  III  ,  page  127 

•1  5AAF  4/18  Himmelblau  ( 1970)  ;  model  8.2-3  ,  page  183 

5  6AAA  3/13  Himmelblau  (1970)  ;  model  6.2-4  ,  page  188 

6  8  AC  A  4/24  None 

7  9AAG  8/25  Inghold  Hertel;  Microwave  Absorption  Line  Function 

(personal  communication) 

8  11AAB  4/9  Ratkowsky  ( 1983)  ;  model  4.4  ,  page  62 

9  12AAB  4/9  Ratkowsky  ( 1983)  ;  model  4.14  ,  page  77 

10  14ACG  3/LO  Ratkowsky  (1983)  ;  model  3.5  ,  page  51  and  58 

11  11ABG  3/21  Rat  kowsky  (1983)  ;  model  3.5  ,  page  51  and  58 

12  14AAG  3/42  Ratkowsky  (1983)  ;  model  3.5  ,  page  51  and  58 

13  15AAA  3/16  Ratkowsky  (1983)  ;  model  6.11  ,  page  120  and  58 

14  16AAF  5/27  Ratkowsky  (1983)  ;  model  6.12  ,  page  122,  123  and  125 

15  17 AAA  2/42  Ratkowsky  ( 1983)  :  model  3.4  ,  page  50  and  58 

16  18AAA  3/9  Ratkowsky  (1983)  ;  model  4.1  ,  page  61  and  88 

17  19AAA  3/9  Ratkowsky  (1983)  ;  mode!  4.2  ,  page  61  and  88 

18  20AAG  4/9  Ratkowsky  ( 1983)  ;  mode!  4.3  ,  page  82  and  88 

19  21  AAA  4/9  Ratkowsky  ( 1983)  ;  model  4.5  ,  page  63  and  88 

20  22AAB  3/5  Ratkowsky  ( 1983)  ;  model  5.1  ,  page  93  and  102 


1 

method  i>c*-or-rrr  ukclmooo  -  uncadiiatioh  - 

VAJHAHT  A  VABUU.T  •  VAAIAMT  C 


Observed  coverage  for  nominally  95  9o  confidence  regions 
with  e~  A>'(0,ct2  I) 
versus 

method  by  derivative  configuration 
Figure  1 


METHOD:  LAC*-Of-m  UCtLHOOO  - - —  t«t»W»TIOH  - 

V  ARtAKT  A  VAAtANT  B  VAAUNT  C 


Observed  coverage  for  nominally  959o  confidence  regions 
with  £-n(o,(tj  a)- 1) 


versus 

method  by  derivative  configuration 


MC’XOO 


L4CB-or-rrr  u«u*ooo  * -  u*ca»uat>ow  - 

VMUkllT  A  VAAIANT  t  vajhamt  c 


Observed  coverage  for  nominally  95fo  confidence  intervals 
with  e~  A'(0,<t2  I) 
versus 

method  by  derivative  configuration 
Figure  3 


MCTWOO: 


lAOt-Or-HT  UttUMOOO 


— -  UNCAAtfATOW 

VARIANT  A  V AQUMT  ■  VAIIIANT  C 


Observed  coverage  for  nominally  Qo^o  confidence  intervals 
with  c-A'(0,(-n  d):  i) 


versus 

method  by  derivative  configuration 


Figure  4 


SCALED  NTffWSC  CURVATURE 

Likelihood  method  confidence  region  observed  coverage 

versus 

intrinsic  curvature  scaled  by  (Ff  i0  os)~I/2 


Figure  6 


END 

FILMED 

10-85 

DTIC 


